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ABSTRACT 

By solving the time- dependent radiation transfer problem of stellar radiation 
through evolving superbubbles within a smoothly varying H I distribution, 
we estimate the fraction of ionizing photons emitted by OB associations that 
escapes the H I disk of our Galaxy into the halo and intergalactic medium 
(IGM). We consider both coeval star-formation and a Gaussian star-formation 
history with a time spread a t = 2 Myr. We consider both a uniform H I 
distribution and a two-phase (cloud/intercloud) model, with a negligible filling 
factor of hot gas. We find that the shells of the expanding superbubbles quickly 
trap or attenuate the ionizing flux, so that most of the escaping radiation 
escapes shortly after the formation of the superbubble. Superbubbles of 
large associations can blow out of the H I disk and form dynamic chimneys, 
which allow the ionizing radiation to escape the H I disk directly. However, 
blowout occurs when the ionizing photon luminosity has dropped well below the 
association's maximum luminosity. For the coeval star-formation history, the 
total fraction of Lyman Continuum photons that escape both sides of the disk 
in the solar vicinity is (f eS c) ~ 0.15 ± 0.05. For the Gaussian star formation 
history, (/ eS c) ~ 0.06 ± 0.03, a value roughly a factor of two lower than the 
results of pove fc Shull ( 1994b )| , where superbubbles were not considered. For a 
local production rate of ionizing photons ^L y c — 4.95 x 10 7 cm -2 s _1 , the flux 
escaping the disk is $L y c ~ (1-5 — 3.0) x 10 6 cm -2 s _1 for coeval and Gaussian 
star formation, comparable to the flux required to sustain the Reynolds layer. 
Rayleigh- Taylor instabilities exist early in the OB association's evolutionary 
stages, possibly causing the shell to fragment and increasing (/ osc ). However, if 
a significant fraction of H I is distributed in cold clouds with nu ~ 30 cm -3 , 
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(fesc) can be reduced by a factor of ~ 2 — 5 if the cloud properties are similar to 
"Standard Clouds" with a disk geometry. 



Subject headings: H II regions — interstellar medium: diffuse ionized gas - 
radiation transfer: photoionization — 



1. Introduction 



The diffuse, ionized medium (DIM), a.k.a the "Reynolds layer," has been established 
to occupy a significant fraction of the interstellar medium (ISM) and to have a vertical scale 
height of ~ 1 kpc ( jMezger 19781 |Keynolds 1991a| ; |Keynolds 199 1 b| ) . It is widely believed 
that OB associations are the sources of radiation responsible for ionizing the Reynolds 
layer ([Reynolds 1984| ; |Mathis 1986| ; |Bregman fc Harrington 1986|) . It is still unclear how 
this ionizing radiation is able to travel so far from the OB associations' immediate vicinity. 
In a previous paper ( Pove fc Shull 1994b| , hereafter referred to as DS94), we calculated 
the geometry of diffuse H II regions due to OB associations in the Galactic plane by 
assuming a three-component model ( |Dickey fc Lockman 1990]) for the vertical distribution 
of hydrogen. For associations producing an ionizing luminosity Sq ^ 3 x 10 49 s _1 , we found 
that the H II regions are "density bounded" in the vertical direction, allowing photons 
to escape into the halo of our Galaxy. By integrating over the luminosity function of 
associations, dN a (So) / dSo ( [Kennicutt, Edgar fc Hodge 1989| ; |McKee 1997|) , we estimated 
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that roughly 7% of ionizing photons, corresponding to a number flux $L y c = 2 x 10 6 cm 
s _1 , escape each side of the H I disk layer in the solar vicinity and penetrate the DIM layer. 
If this radiation is absorbed in the DIM, assumed to lie entirely above the HI layer, the 
corresponding photoionization rate would be comparable to the recombination rate implied 
by the Ha emission measurements QMezger 1978 ; Reynolds 1991a ; Hcynolds 1991b ). If 
roughly one-third of this radiation escapes the Galactic halo, this flux is also consistent with 
the estimated flux required to photoionize the Magellanic Stream and several high velocity 
clouds in the solar vicinity ( Bland-Hawthorn & Maloney 199"9] ). Rather than assuming a 
purely diffuse ISM, [Miller fc Cox (1993]] considered the effects of a statistical distribution of 
small opaque clouds, also decreasing with increasing Galactic height and embedded within 
the diffuse gas. These authors found that the resulting H II regions are consistent with the 
observations of dispersion measures and Ha emission measurements. 

Although these results are encouraging, no models to date consider the role of dynamic 
chimneys and superbubbles produced by the OB associations. These considerations are 
the focus of this paper. We recently became aware of a preprint on a similar topic (|Basu 
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Johnstone, & Martin 1999). For an OB association, the stars that produce most of the 
ionizing radiation have strong stellar winds and form supernovae after their relatively short 
lifetime. Therefore, shortly after the formation of an OB association, a superbubble will 
form, with a thin shell behind the shock front excavating a large cavity of hot gas (|Weaver 



let al. 19771 ; McCray fc Snow 19791 ; McCray fc Kafatos 198ft IShull fc Saken 1995Q . For 



a sufficiently large association, the superbubble can break out of the H I disk, forming a 
dynamic chimney, where the hot cavity penetrates into the halo (|Mac Low fc McCray 1988| ) . 

As we discuss below, the shells of the expanding superbubbles quickly trap the ionizing 
radiation, so that no radiation can escape the disk until blowout. After blowout, a large 
fraction (~ 10% per side) of ionizing radiation can directly escape the disk of the galaxy 
by propagating through the dynamic chimney. However, this is a viable mechanism for 
transporting the ionizing radiation only if an appreciable amount of ionizing radiation is 
still being produced within the cavity after the dynamic chimney has formed. We find that, 
for both a coeval star formation history or a non-coeval model having a star formation 
rate given by a Gaussian with a full-width half-maximum a t = 2 Myr ( Massey 1998| ; 



Garmany 1998|) , a small fraction of the ionizing radiation produced by an association over 



its lifetime is emitted after the time of blowout. Therefore, most of the radiation that 
escapes the disk does so early in the early stages of the superbubble. Integrating over the 
luminosity function, we find that the fraction of ionizing photons currently being emitted 
by OB associations that escape each side of the H I disk is (/ eS c) ~ 15% for the coeval 
star- formation history and (/ e sc) ~ 6% for the Gaussian star-formation history. 

The aim of this paper is to predict the fraction of ionizing radiation, emitted by OB 
associations, that escapes the H I disk of the Galaxy. As discussed above, DS94 assumed 
that the Reynolds layer lies above the H I layer. Therefore, the question of whether the 
Reynolds layer is sustained via hot star photoionization is equivalent to whether enough 
radiation can penetrate through the H I layer. In this paper, we regard this assumption 
as too simplistic. The Reynolds layer has a scale height of roughly 1 kpc, whereas the H I 
disk has a scale height of ~ 200 pc ( Dickey fc Lockman 1990|) but extends to at least 1 kpc 



( [Lockman, Hobbs, fc Shull 1986| ). Thus, a significant fraction (e.g., ~ 20%) of the of the 
ionized gas may lie within the H I layer, reducing the flux of ionizing radiation required to 
escape the H I disk. A self-consistent model, modeling the distribution of the photoionized 
gas, is not feasible until the density distribution of the hydrogen gas, both within the disk 
and in the halo, is determined more accurately. 

In §2, we discuss our treatment of the evolution of the superbubble and the radiation 
transfer of ionizing photons emitted by an association through the expanding superbubble 
and the diffuse ISM of the disk. We also discuss our treatment of the blowout event of the 
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superbubble, where the cavity evolves into a dynamic chimney. In §3, we determine the 
time dependence of the fraction of ionizing photons, emitted by a single OB association, 
that escape the H I disk. Integrating over the luminosity function of OB associations, we 
determine the fraction of ionizing radiation, currently being emitted by OB associations, 
that escapes the H I disk. In §4, we discuss key issues not included in our standard 
model, including shell instabilities, a two-phase ISM, "poisoning" by photoablated gas, and 
triggered star formation. We then discuss how these effects may alter our results. In §5, we 
discuss the implications of our results and give our conclusions. 



2. Radiation Transfer Through an Expanding Superbubble 

2.1. Evolution of the Superbubble 

Consider a large association having iVr stars in the 8 — 85 M mass range. Due 
to the stellar winds from the massive stars and type II supernovae, the association 
produces a mechanical luminosity, L mec h(t), that varies with time. This energy input 
drives an expanding superbubble. Assuming an initially uniform density distribution of the 
diffuse ISM, and that all of the swept-up mass resides in a thin shell ( the Kompaneets 
approximation, [Kompaneets 19601) , the radius, i? s h, of the cavity is given by flShull fc Saken 
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L m ech(t), (1) 

where the mass density p = lAm p nH, uh is the number density of hydrogen, and m p is the 
proton mass. For a given time-dependent mechanical luminosity prescription, this equation 
is integrated numerically using a fourth-order Runge-Kutta algorithm. We note that this 
formalism is technically valid only for a uniform density distribution, as the shock wave 
does not evolve with spherical symmetry for non-uniform density distributions. However, as 
we show in §3, for most association sizes the expanding shell becomes optically thick to the 
ionizing radiation before the shock has propagated a distance of a disk scale height, within 
which using a constant density is a good approximation. Therefore, the amount of escaping 
ionizing radiation is insensitive to the errors produced by neglecting the non-spherical 
propagation that occurs at large Galactic heights. 

The time dependence of the mechanical luminosity depends on the initial mass function 
(IMF) of stars within an OB association and the star formation history of the association. 
For this paper, we use the OB association evolution models of [Sutherland fc Shull (1999)] . 
which determine both the mechanical luminosity and the luminosity of ionizing photons as 
a function of time for individual stars from the latest stellar evolutionary tracks and stellar 
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atmospheres. These authors used an IMF given by dN/dM oc M~^ +l \ where 7 = 1.6 for 
the stellar mass range 8 M < M < 85 M ( parmany 1998| ), for which all stars eventually 



become supernovae. The massive star-formation history of the association is given by 
N*(t), with the constraint that J iV*(i)<it = Nt, where the limits of integration span the 
entire lifetime of the association. In this paper, we consider two star-formation histories 
(SFHs): (1) coeval [N*(t) oc 5(t — t )], and (2) noncoeval with a Gaussian distribution 
[N*(t) oc exp(— £ 2 /2of )]. For the Gaussian SFH, the IMF is assumed to be constant in time. 



For each SFH, L mec h(£) and S(t) were obtained from [Sutherland fc Shull (1999 



Here, we used a Monte-Carlo simulation to determine the mass distribution of individual 
associations having a total mass of 10 4 M (corresponding to an average association size of 
3000 stars, or roughly ~ 300 massive stars). The average luminosity curves are then 
normalized, and, for an arbitrary association size, the luminosity curves are given by the 
product of the normalized curves and iVr- Note that this formalism ignores all possible 
deviations about the average luminosities due to statistical fluctuations of the stellar mass 
distribution of small associations. Sutherland & Shull (1999)1 found that these fluctuations 



become significant for iVr <^ 100. However, as shown below, most of the escaping radiation 
in our Galaxy emanates from associations with ^ 100, so the importance of the 
statistical fluctuations of small associations is minor for our purposes. Even if the small 
associations were important, when calculating the total fraction of escaping photons, we 
integrate over many associations, so the use of the average luminosities should yield a good 
approximation. 

In Figure |], we show the mechanical luminosity and the Mach number of the shell as a 
function of time for an OB association having = 200, for coeval star-formation history 
and for a Gaussian star-formation history with o t = 2 Myr. The Mach number of the shock 
front is M. = f s h/c s , where the effective sound speed of the medium is 



/,T -(l+/3)- (2) 



// p m p 

Here, k is the Boltzmann constant, T is the temperature of the gas (assumed equal to 10 4 
K and the shock front is assumed to be isothermal), fi p = 0.6, and /3 is the ratio of the 
magnetic pressure plus turbulent pressure to the thermal pressure. As discussed in § pT2| , the 
superbubble is modeled as instantaneously blowing out of the disk and becoming a dynamic 
chimney when the bubble radius exceeds the transition height, Z tr . We assume that Z tT is 
twice the scale height of the H I gas distribution, a h . After blowout, the cylindrical radius 
continues to expand due to conservation of momentum. Once the Mach number reaches 
unity, the expansion of the cavity is halted. 

The fluctuations in the mechanical luminosity arise from supernovae occurring 
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Fig. 1. — Mechanical luminosity of the OB association and Mach number of the shell as a 
function of time. For the Gaussian SFH, a t = 2 Myr and the star-formation peak occurs 
at t — 2<jf For both cases, iVr = 200 and (3 = 1. Blowout occurs when t ~ 11 Myr for 
the coeval SFH and t ~ 16 Myr for the Gaussian SFH. After blowout, the shell is in the 
momentum-conserving phase, and the Mach number shown is evaluated at the midplane. 
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separately in time due to the different lifetimes of individual stars. Although L mec h(£) is 
given by the average of several Monte-Carlo simulations, the fluctuations have not been 
smoothed out. 

We note that, for the Gaussian SFH, the Mach number is near zero until 1.3 Myr, 
when it then rises quickly. This behavior is due to our initial treatment of the bubble 
growth for low mechanical luminosities. Since we start the simulation at a time 2a t before 
the peak star-formation rate, the initial luminosity is small for small associations. In such 
cases, the Mach number rapidly drops below unity, then subsequently increases as the 
energy input increases. In reality, the shock front associated with these small associations 
would stagnate when the ambient pressure of the ISM is equal to the internal pressure 
of the bubble, causing the Kompaneets approximation to break down. As the proper 
hydrodynamic treatment of the growth of a stagnated bubble is outside the scope of this 
paper, we employed the following procedure for approximating the bubble evolution. For 
each association size, we determine the time at which the mechanical luminosity exceeded 
10 36 erg s -1 . We then evaluate the Mach number at this time using the solution of equation 
(0) for a constant luminosity, where the time-averaged luminosity is used. If the Mach 
number is below unity, we assume that the shock has stagnated, set the shell velocity to 
zero, and continue integrating equation ([]]). Although the use of Kompaneet's equation is 
not valid for Mach numbers below unity, we argue that the errors introduced are negligible 
since the time interval for the shell to accelerate to Mach numbers above unity is small 
relative to the sound-crossing time of the bubble. Therefore, the shell does not diffuse much 
into the cavity during the stagnation period. Furthermore, during this early interval of the 
cavity, the photon luminosity of the association is also low, and no radiation escapes the 
disk of the galaxy regardless of the bubble geometry. Finally, by comparing M. (t) using this 
method with results when equation (|l|) is integrated without intervention, we found that 
the two were identical for t ^ 1.5 Myr. The only ramification of our stagnation treatment 
is that smaller associations take longer to blow out of the disk. As discussed in §fO, the 



acceleration of stagnated bubbles due to the mechanical luminosity increasing in time can 
cause Rayleigh- Taylor instabilities to begin at a much earlier stage, compared to the case of 
constant luminosity, where acceleration occurs at the time of blowout ( Mac Low fc McCray 



1988|) . This instability may dramatically increase the amount of ionizing radiation that 



escapes the H I disk. 
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2.2. Superbubble Blowout 

For large OB associations that have a sufficiently high mechanical luminosity, 
superbubbles can blow out of a vertically stratified H I disk, producing dynamic chimneys 
( gales 1979 1 ; ITomisaka fc Ikeuchi 1981 ; |Mac Low fc McCray 198$ |Mac Low, McCray & 



|Norman 1989Q . Blowout occurs if the shell of the superbubble is still supersonic when it 
reaches a height, Z acc , where it begins to accelerate due to decreasing density distribution 
( [Mac Low fc McCray 1988|) . This height depends on the density distribution of the ISM and 
L mec h', for Gaussian or exponential distributions it is Z acc ~ (2 — 3)o"h (|Mac Low fc McCray 



1988fe |Mac Low fc Ferrara 1999| ). For most association sizes and star formation histories, the 



timescale for reaching Z acc is considerably longer than the timescale for propagating from 
this height into the Galactic halo [Z ^ 1 kpc). Therefore, rather than calculating the exact 
evolution of the superbubble during blowout, we model the process as an instantaneous 
transition from a bubble into a dynamic chimney once the shell reaches a transition height 
Ztt = Z acc . We note that, due to the existence of the H I gas at high Z, blowout may 
first occur at heights considerably larger than 3<7h or may not occur at all. In addition, 



magnetic fields could inhibit blowout flTomisaka 1998|) . As discussed below, our results are 



insensitive to the exact value of Z tr , since most ionizing radiation that escapes the disk does 
so shortly after the formation of the OB association. Therefore, the uncertainties regarding 
the process of blowout are unimportant for our estimation of the flux of escaping ionizing 
radiation. 

The geometry of the dynamic chimney is approximated as a cylinder, centered on 
the OB association, with radius R cy and height Z CJ . Since we are only interested in how 
much ionizing radiation escapes the H I disk of the Galaxy, we do not attempt to model 
the structure of the superb ubble/chimney for heights more than ~ 4 — 5 disk scaleheights 
1 kpc). Numerical simulations of superbubbles blowing out of stratified disks ( [lbmisaka 



[I998t |Mac Low, McCray fc JNorman 1989j ; |Mineshige, Shibata fc Shapiro 1993|) show that a 



cylinder is a decent approximation for the geometry of the cavity within the H I disk. Above 
the disk, the shell takes the appearance of a mushroom cloud, with the shell clumping 
and fragmenting due to Rayleigh- Taylor instabilities. For this paper, we assume that any 
radiation that reaches a height Z cy and is within the dynamic chimney escapes the disk. 
The question of whether this radiation is absorbed by the "mushroom cap" of the shell or 
escapes due to the shell having a small covering fraction is outside the scope of this paper. 
Even if the covering fraction of the "mushroom cap" is not negligible, it is possible that 
this gas, photoionized by the OB association, is part of the Reynolds layer. Therefore, it is 
fair to include this radiation when assessing whether the amount of escaping radiation can 
ionize the DIM. 
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2.3. Ionization Equilibrium of the H II Region 



Throughout this paper, we assume that ionization equilibrium is satisfied 
instantaneously. As discussed in DS94, the ionization front in a Gaussian disk layer 
will reach a Galactic height Z in a time 



where 649 is the photon luminosity is units of 10 49 s -1 , is the scaleheight of Gaussian 
disk, and no is the number density of hydrogen in the midplane of the disk. (This equation 
should replace equation (3) of DS94, whose equation (3) contains an error.) For S49 = 1, 
o"h = 0.184 kpc, and n = 0.367 cm -3 , the I-front reaches a height Z = in a time t ~ 1.4 
Myr and reaches a height Z = 2<7h in a time t ~ 3.2 Myr, both relatively small compared to 
the lifetime of the superbubble. As shown in DS94, ionization equilibrium is reached within 
a decade or so after the passing of the I-front. [For gas with an ionization fraction X ~ 1, 
ionization equilibrium is reached in a time r ~ (1 — X) / {n e a^)\. Thus, even though the 
density distribution of the superbubble is evolving in time (as discussed below), the H II 
region is determined in a quasi-static fashion. 



We assume that the stars of individual OB associations are compactly distributed so 
that the radiation can be treated as emanating from a point source. We define y(t) to 
be the dimming function of the OB association, which relates the time-dependent photon 
luminosity of Lyman-continuum (LyC) photons and the maximum luminosity via 



For a point source of ionizing radiation, the H II/H I boundary is given implicitly by 



t(Z) 




(3) 



2.3.1. Geometry of the H II Regions 



S(N T ,t)=S {N T )y(t). 



(4) 








(5) 
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where 9 is the angle between the normal of the Galactic disk and the radial unit vector, a H 
is the case-B recombination rate for hydrogen, and n# is the number density of hydrogen. 
Azimuthal symmetry has been assumed. 

For the vertical distribution of H I, DS94 considered both the three- component 
Dickey-Lockman model ([Dickey fc Lockman 1990|) as well as several single Gaussian 
distributions (with different amounts of dark matter) which were fitted to numerical 
calculations assuming hydrostatic equilibrium ( Pove fc Shull 1994a|) . We found that 
all models predicted the same fraction of escaping photons within 10%. This Gaussian 
distribution which most closely resembled the Dickey-Lockman model is given by 



iih{Z) = n exp(— Z 2 /2a i 



(6) 



where Z is the vertical distance above the mid-plane of the disk, n = 0.367 cm 3 and 
ah = 0.184 kpc. In this paper, we only consider this single- Gaussian distribution. 

As the superbubble expands into the diffuse ISM, all of the swept-up mass is assumed 
to reside in the thin shell. The thickness of the shell, AR, is given by equating column 
densities, 

^sh 

n sh (R sh ,fx)AR= J n H {Z)dR, (7) 

Q 

which yields 



AR=*ll(?!L) erf (**L) 
n sh \2\fxj \V2aJ 



where n s h is the number density within the shell, assumed to be a constant with respect 
to radius, and fi = cos 8. Defining the transition radius R tT = i? s h + A_R, the geometry of 
the H II region for an OB association with luminosity S, situated on the mid-plane of the 
Galaxy, is given implicitly by 



S 



Ana 



(2) 
H 



- t nl h (R sh , fi)(Rl T - R, 



l sh) 



(9) 
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The density of this shell is given by the isothermal shock condition for compression in the 
cool post-shock layer, 

n sh {R sh , fi) = M 2 n R {fiR sh ). (10) 

In Figures [| and 0, we show the time evolution of the geometry of the H II region 
shortly after the formation of an OB association, for = 50 and = 300, respectively, 
assuming coeval star formation. For both association sizes, the H II region is initially 
density bounded in the vertical direction (see below for discussion of how much radiation 
escapes the H I disk). As the superbubble expands, the column density of the shell 
increases. Since the shell has a higher density than the diffuse gas, and therefore a higher 
recombination rate, the volume of the H II region decreases with time. Eventually, the 
H II region is radiation bounded even in the vertical direction. For = 300, a portion of 
the H I/H II interface lies within the shell for t ^ 3.0 Myr, and by t ~ 3.6 Myr the entire 
interface is within the shell. The H II region goes from being density limited to radiation 
limited in a short time interval. For = 50, the H II region becomes radiation limited at 
only a slightly earlier time compared to the = 300 association (2.9 Myr compared to 
3.5 Myr) even though its photon luminosity is six times lower. This is due to the smaller 
association having lower shell velocities, causing both a slower accumulation rate of the 
column density as well as lower gas densities within the compressed shell. 



2.3.2. Escaping Radiation 



If the H II region is density bounded, then some radiation escapes the H I disk. The 
fraction of ionizing photons that escape along the differential solid angle dQ(9) = sin# d9dc l 
is 



f[e,s(N T ,t)\ 



Ana 



(2) 
H 



S(N T ,t) 



n 2 H (R,9,t)R 2 dR 



S(AM)1 3 Uo ' 



<7h 
fi 

2a h 



1 — erf 



+ 



■ exp 



tr 



yR 

Oh 



(11) 



The critical angle, C (S), is the angle at which the H II region is radiation bounded for 
9 > 6 C and density bounded for 9 < 9 C . The critical angle is found by solving f(9 c , S) = 0. 
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Fig. 2. — Geometry of H II region for several times after the formation of a coeval OB 
association with = 50 massive stars. From right to left, the time (in Myr) after formation 
is t = 0.0, 1.6, 2.2, 2.4, 2.6, 2.7, 2.8, 2.85, 2.9, and 4.0. 
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Fig. 3. — Geometry of H II region for several times after the formation of a coeval OB 
association of size iVp = 300. From right to left, the time (in Myr) after formation is 
t = 0.0, 2.20, 2.77, 2.93, 3.13, 3.36, 3.40, 3.50, 3.51, 3.54, and 3.61. 
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In contrast to DS94, 9 C is time-dependent since both the density distribution and the 
ionizing luminosity vary with time. Finally, we define the fraction of photons, emitted by 
a single OB association having size iV T , age t, and corresponding luminosity S(N T ,t), that 
escape each side of the H I disk to be 

9c 

Vesc(N Tl t) = I dO. (12) 



2.3.3. H II Regions for Dynamic Chimneys 

For a dynamic chimney with height Z CJ and cylindrical radius R cy , the angle for which 
f{9) = 1 for < 9 < 9d yn is given by 9d yn = tan -1 (i2 cy /Z cy ). The procedure for determining 
the geometry of the H II region and the fraction of escaping photons is similar to that of the 
superbubbles. The only difference is the treatment of the shell. As with the superbubble, 
we assume that the density of the shell is n s h = M 2 nn and that the column density of the 
shell is equal to the column density of swept-up gas within the cavity. However, for dynamic 
chimneys, the trajectories of the shell deviate from radial lines since the superbubble 
evolves from a sphere into a cylindrical cavity. Technically, the column density of the shell 
for a vertical height Z would be given by the line integral of the density of the diffuse ISM 
along the trajectory that leads to the chimney wall at a height Z. The trajectories can be 
determined by requiring them to be perpendicular to the surface of the expanding shock 
wave. For an isothermal atmosphere, [Newman et al. (1999)| found an analytic solution for 



the trajectories as a function of time, and, for heights less than a scale height or so, the 
particle trajectories are found to be roughly radial. For larger heights, the trajectories 
tend to bend away from the vertical axis and become horizontal near the shell boundary. 
Therefore, we make the approximation that the streamlines are purely radial for Z < Z tr 
and purely horizontal for Z > Z tI , where Z tr ~ 2 — 3(Xh, a free parameter of the model. 
As discussed below, the results of this paper are insensitive to the value of Z tT . With this 
treatment, for Z < Z tr , the shell thickness is given by 

i 

A R( 7\ f 

n s h — = / n n dR, (13) 



sin 9 

b 

where I = ^K 2 cy + Z' 2 and 9 = tan' 1 (R cy / Z) . Defining R sh = I and R tI = I + AR(Z)/sin9, 
the geometry of the H II region and the fraction of escaping radiation are again given by 
equations (P) and (0), respectively. 
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3. Results 

3.1. Fraction of Escaping Photons as a Function of Time 

3.1.1. Coeval Star Formation 

In Figure f| we show the photon luminosity, emitted by a single OB association having 
iVr stars, and the luminosity of photons escaping each side of the H I disk as a function of 
time. Here, the stars are assumed to form coevally. Notice that the fraction of escaping 
photons, 77esc(^Vr, £), rapidly decreases for t ^ 3 Myr. In fact, for N T = 1000, roughly 90% 
of the time- integrated flux of escaping radiation occurs within the first 2.7 Myr, while only 
72% of the time-integrated flux of ionizing radiation has been emitted by this time. The 
reason for this rapid decrease is that, as the column density of the shell accumulates, the 
ionizing radiation is very quickly trapped within the shell, causing the fraction of escaping 
photons to drop to zero. Once the ionization front is trapped within the shell, the fraction 
of escaping photons remains zero unless the superbubble blows out of the H I disk. (We 
have not considered shell fragmentation). If blowout occurs, photons can directly escape 
through the chimney (a small fraction can also escape by penetrating through the upper 
H II disk), causing r] esc to jump up to ~ 10%. However, by the time of blowout, which 
occurs at t ~ 7.0 Myr for = 1000 and 10.0 Myr for = 200, the photon luminosity 
is considerably lower than its peak luminosity (which occurs at t ~ 2 Myr). Therefore, the 
luminosity of escaping radiation is significantly lower compared to its value for t <^2 Myr. 

3.1.2. Gaussian Star Formation 

In Figure |5| we plot the emitted and escaping photon luminosity luminosity for case of 
the noncoeval, Gaussian SFH with cr t = 2 Myr. For t = 0, the photon luminosity S(t) is 
very low, and therefore the fraction of escaping photons is small or even zero, in accordance 
with the results of our DS94. As the luminosity increases, r] esc increases until the shell of the 
superbubble becomes thick enough that the shell-integrated recombination rate approaches 
the photoionization rate. Since the shell expands more slowly compared to the coeval SFH 
(having the same value of N^), the column density of the shell does not increase as quickly. 
In addition, due to the lower Mach number of the shock front, the shell density is lower 
than that of the coeval SFH. For these two reasons, i] esc does not fall off in time as quickly 
as in the coeval case after it reaches its peak value. Nevertheless, the shell does attenuate 
significantly the flux of escaping radiation. In Figure for = 200, we compare the flux 
of escaping radiation as a function of time with that for the case where no dynamics are 
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Fig. 4. — Photon luminosity emitted by a single OB association (solid line) and the 
luminosity of photons escaping each side of the H I disk (dashed line) for coeval star- 
formation. Also shown is the fraction of photons emitted that escape, ^^(N^jt) = 
S eS c(t)/S(t). Here, Z tr = 2a^, the cavity height is Z cy = 4<7h, and (3 = 1. The I-front 
is trapped within the shell during the time interval where ?7esc = 0. After blowout, ?7esc 
suddenly increases even though the photon luminosity of the association is relatively low at 
these times. 
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considered (e.g., as done by DS94 but using the time-dependent photon luminosity). For 
t ^ 5 Myr, the shell becomes becomes efficient in reducing the escaping flux. Therefore, 
as with the coeval SFH, most of the escaping radiation escapes early, within a short time 
interval relative to the lifetime of the OB association. The time-integrated flux of photons 
that escape the disk after blowout is considerably small relative to the lifetime-integrated 
flux. For example, for = 1000, fewer than 3% of the photons that escape the disk over 
the lifetime of the association escape after blowout (while fewer than 4% of the photons 
emitted by the association are produced after blowout). For iVr = 200, approximately 1% 
escape after blowout. These values are not much different than for the coeval SFH. The 
reason for this insensitivity is that as a t is increased (er t = for coeval star- formation), the 
blowout timescale is also increased. Thus, even though there are more stars alive at longer 
times for noncoeval SFHs, there are still relatively few alive at the time of blowout. 



3.2. Time Integrated Fraction of Escaping Photons 

The time-integrated fraction of escaping photons for each side of the H I disk by an 
association of size is given by 

oo 

/ r] esc (N T ,t)S(N T ,t) dt 

(vUNt)) = ^ - . (14) 

J S{N T ,t) dt 
o 

In Figure |[] we plot (t^sc^t)) for both the coeval SFH and the noncoeval Gaussian 
star-formation model. For reference, the maximum photon luminosity of an association over 
its lifetime is S (N r ) = 0.231iV T (10 49 s" 1 ) for coeval star- format ion and S (N T ) = 0.122iV T 
(10 49 s™ 1 ) for the Gaussian SFH. Note that (r] esc ) is smaller for all values of iVr for the 
Gaussian SFH. There are several reasons for this difference. One reason is that, for a given 
value of iVr, the Gaussian SFH has a peak photon luminosity that is roughly a factor of 
two smaller than that of the coeval model. Therefore, even if the shell structures of the 
superbubble were identical during the interval of maximum luminosity, the maximum value 
of 7/esc for the Gaussian model would be smaller since r] esc increases with increasing photon 
luminosity (DS94). In addition, the luminosity profile S(t) for the coeval model is much 
more peaked, with ~ 70% of its lifetime-integrated luminosity emitted within the first 2.5 
Myr, compared to the relatively broad distribution of S(t) for the Gaussian SFH. Thus, 
since (r} esc ) is weighted by S(t), and t^sc^t) is near its peak value while S(t) is near its 
peak value, the coeval SFHs yield higher values. Finally, for coeval SFHs, the peak of S(t) 
occurs shortly after the formation of the association, before the formation of a significant 
superbubble shell, allowing these photons to escape more easily. 
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Fig. 5. — Photon luminosity emitted by a single OB association (solid line) and the 
luminosity of photons escaping each side of the H I disk (dashed line) for a Gaussian 
SFH with at — 2 Myr. Also shown is the fraction of emitted photons that escape, 
Vcsc(N T ,t) = S esc (t)/S(t). Here, Z tT = 2a h , Z cy = Aa h , and (3 = 1. 
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Fig. 6. — Photon luminosity emitted by a single OB association (solid line) of size = 200 
and the luminosity of photons escaping each side of the H I disk (dotted line) for a Gaussian 
SFH with a t — 2 Myr. Also shown is the luminosity of escaping photons for the case where 
the dynamics of the bubble is not considered (dashed line). The bottom panel shows the 
Mach number of the shock front as a function of time. Due to the high density of the evolving 
shell, the luminosity of escaping radiation is reduced for t ^ 4 Myr. 
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Fig. 7. — Time-integrated fraction of escaping photons as a function of association size iVp. 
Here, Z tT = 2<r h , Z cy = 4<r h , and (5 — 1. Solid line: Gaussian SFH, a t — 2 Myr. Dashed line: 
Coeval SFH. 



3.3. Fraction of Escaping Ionizing Photons 



We now discuss the fraction of photons, currently being emitted by OB associations of 
all ages and sizes, that escape the H I disk of the Galaxy. 



3.3.1. Luminosity Function 

We assume a constant formation rate, iV as , of associations per unit area in our 
Galaxy. The probability of a new association having a peak luminosity Sq [recall that 
S(NT,t) = So(Nj:)y(i)] is given by the implicit association luminosity function, dN/dSo- 
Therefore, during a small time interval 8t, the number of associations that will eventually 
have a peak luminosity S is 5N(S ) = N^dN/dSo St. The number density (kpc~ 2 s _1 ) 
of associations observed today as a function of luminosity, dN as (S)/dS, is related to the 
implicit luminosity function via 

dN as (S) _ f . v dN 



ijs j AT as — [So = S/y{t)\ dt. (15) 

h(s) 

Here, t 2 (S) is given by t 2 = y~ 1 (S/S 2 ), where y~ l is the inverse function of y(t), and t\ is 
given by 

* -! : S>S 1 

^-{y-^S/S,) : S<S l . [ib) 

If the dimming function y(t) is multi- valued, then t 2 and t\ (for S < S\) are determined 
from its maximum values. 

The observed luminosity function, from Ha measurements of all known H II regions in 
30 nearby spiral and irregular galaxies, is given by dN as {S)/dS oc S~ r for S\ < S < S 2 , 
where T = 2.0 ± 0.5 ( [Kennicutt, Edgar fc Hodge 1989| ; Banfi et al. 1993| ; |Rozas, Beckman"^ 



Knapen 1996| , and references therein). [McKee (1997)| finds a similar result from studying 



Galactic H II regions. See DS94 for a discussion regarding the uncertainty of the upper 
and lower limits of the luminosity function. Assuming dN/dSo oc Sq F with T ^ 1, the 
corresponding dN as (S)/dS is close to a pure power law with the same power-law index, 
with deviations first occurring at S ~ S 2 . Therefore, for the remainder of this paper, unless 
otherwise stated, we assume that the intrinsic luminosity function is a power law with 
T = 2. 



A better understanding of the relationship between the implicit luminosity function 
and the observed luminosity function can be obtained by considering an example with 
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a simple analytical dimming function, y(t). Assume y(t) = exp(— t/r) and the implicit 
luminosity function dN/dS = N Sq T . Then, for S > S 1: 

t 2 (S) 

dN as (S) f aV dN 



dS 



J K^(S = S/y(t))dt (17) 



o 



-r 



= N^NoS' 1 ' J exp(-n/r) dt 
o 

oc ,S- r [l-exp(-rt 2 /r)]. (18) 
Here, t 2 = y~ 1 (S/S2) = —t\yl(S/ S2), so finally we have 

dN^(S) 



dS 



oc S~ T 



1 - 



5 2 



(19) 



Therefore for S > S\ and S <^ S2, dN^/dS oc S~ r . The larger the value of T, the closer 
S must be to S2 for appreciable differences to occur between the intrinsic and observed 
luminosity functions. For T = 2, a 10% difference between the two functions first occurs 
when S/S 2 > 0.32. Since the range of association luminosities is several orders of magnitude, 
and the value of S 2 is very uncertain, we regard this discrepancy as negligible. Note that the 
reason for any differences between the two luminosity functions is due to our prescription 
that the implicit luminosity function is truncated at S 2 - Physically, as S approaches S 2 , 
only the largest associations recently formed can contribute to the luminosity function, 
causing the observed luminosity function to drop to zero as S reaches S 2 . Also note that 
the difference between the two luminosity functions is independent of the functional form 
of the dimming function. 



3.3.2. Production Rate and Escaping Rate of Ionizing Photons 

The current production rate of ionizing photons (# per unit area per unit time) is 
given by 

* LyC = / s ^nr dS - (20) 
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Substituting the expression for the observed luminosity function in terms of the implicit 
luminosity function [equation ([T5D1, we find, 



5 2 H 



^ Ly c = ^f S J Jj^ S * = S/y(t)) dt dS. (21) 



t\ 

The flux of photons escaping from each side of the H I disk is given by 

5 2 H 



/p 1 71 T 

S I Vosc(S,t)— [S = S/y(t)} dt dS. (22) 



o t 



Here, r) eBC (S,t) is the fraction of photons that escape each side of the H I disk 
emitted by an association of size iVp, where the peak luminosity of this association is 
max[S'(A r T , t)] = S (N T ). By computing a two-dimensional grid of r] esc (N T ,t) using the 
methods outlined in § [2.3.2| , we integrate equation (^) for the two star-formation scenarios. 
For given values of S and t, an interpolation routine was used to find the corresponding 
value of iV T , such that rj esc (S,t) = r] esc [N T (S,t),t}. Finally, the total fraction of photons, 
currently being emitted, that escape both sides of the H I disk is 

</e S c) = 2|=L. (23) 

We find that (/ csc ) ~ 12% for the coeval SFH and (/ esc ) ~ 6% for the Gaussian SFH. 
In Figure we show how (/esc) depends on the model parameters. For reasons given in 
§3.2j , (/esc) for the Gaussian SFH is considerably lower than that for the coeval SFH. By 
far, (/esc) is most sensitive to the power-law index of the luminosity function, T. This 
sensitivity arises because most of the escaping radiation comes from the relatively few large 
associations with iVx ^ 100; small changes in T cause rather large changes in the relative 
number of these large associations. For the coeval SFH, (/ esc ) is insensitive to the magnetic 
field (the parameter (3). As the superbubble expands, the column density of the shell 
increases quickly, causing rj esc (S,t) to drop to zero very quickly even for strong magnetic 
fields. For the Gaussian SFH, the superbubble expands more slowly, and different values of 
(3 cause the shell to trap the ionization front at different times. Since this trapping of the 
I-front occurs near the peak of S(t), (/ csc ) is more sensitive to (3. Lastly, we do not plot how 
(/esc) depends on the parameters involved with the superbubble blowout (e.g., Z tT and Z cy ), 
as we found that varying these parameter over a range between ah and 4a^ corresponded to 
a relative change of (/ esc ) less than 1%. This insensitivity occurs because most emitted and 
escaping photons occur early after the formation of the association. 
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Fig. 8. — Fraction of Lyc photons, currently being emitted, that escape each side of the 
H I disk. Solid line: coeval SFH. Dashed line: Gaussian SFH with a t — 2 Myr. Top panel: 
r = 2 (dN/dSo oc Sq T ) and the maximum association size is N T = 3200 (corresponding to 
^2 = 7.4 x 10 51 s _1 for the coeval star formation history and S 2 = 3.9 x 10 51 s _1 for the 
Gaussian SFH). Middle panel: (3 = 1, and the range of association sizes is the same as given 
above. Bottom panel: T = 2 and (3 = 1. For all models, Z tT = 2a h . 
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4. Complexities to the Model 

We now discuss the key issues not included in our standard model. We then discuss 
how these effects may alter our results. 

4.1. Shell Instabilities 

As discussed above, the presence of a thick neutral shell of cool, shocked ISM produced 
by multiple supernova explosions is quite effective in absorbing the radiation coming from 
the OB association. This process significantly reduces the fraction of photons able to escape 
into the galactic halo until blowout takes place. At that stage, the previously decelerating 
shell is suddenly accelerated and the interface between the dense shell and the hot cavity 
gas becomes Rayleigh- Taylor (R-T) unstable. The nonlinear development of such instability 
induces fragmentation of the shell into a number of clumps, through which the photons can 
percolate and almost freely escape into the halo. This process generally occurs relatively 
late (or not at all), when the radiative flux has already decreased by a considerable amount 
due to the aging of the parent OB association. It is thus worthwhile to investigate if the 
fragmentation process driven by R-T or other instabilities can take place before blowout, 
consequently increasing (/ esc ). 

There are at least three types of instabilities that could affect the shell at the beginning 
or during the pressure-driven radiative (PDR) phase: 1) dynamical, 2) Rayleigh- Taylor 
and 3) gravitational. To assess the importance of these processes, it is necessary to discuss 
in more detail the development of the shell structure for SN bubbles. We restrict our 
discussion to uniform ambient gas density, which is a reasonable approximation to study 
the pre-blowout phase. The structure of a wind-blown bubble (the energy injection by 
multiple SNe can be well approximated as a continuous process) is constituted by a reverse 
( "wind" ) shock propagating back into the wind, which is under most conditions adiabatic, 
and an "ambient" shock propagating into the surrounding medium; the two shocks are 
separated by a contact discontinuity. Initially, the ambient shock is also adiabatic and its 
radius increases with time as t 3 ^ 5 . However, when the cooling time of the shocked ambient 
gas becomes of the order of the age of the bubble, the ambient shock becomes radiative and 
the shocked ambient gas collapses into a thin shell. 

The shell-formation phase is characterized by violent dynamical instabilities which 
tend to amplify the perturbation due to classical thermal instability in the catastrophically 
cooling shocked ambient gas. Are these instabilities able to disrupt the shell soon after its 
formation? At least three different dynamical instabilities may arise in the thin, cool shell: 
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nonradial oscillations of the radiative shock ( [Bcrtschingcr 1956] ), the nonlinear thin-shell 
(NTSI) instability arising when the shell is bounded by two radiative shocks ( |Vishniac 



, and the pressure-driven thin-shell overstability (PDTSO) found in a shell bounded 



by a high thermal pressure on one side and by a radiative shock on the other ( |Vishniac 



|1994|) . These dynamical instabilities may produce considerable distortions in the shell, often 
showing nonradial oscillations. The main difference between the latter two instabilities is 
that while the NTSI grows exponentially, the PDTSO only grows with a power-law rate. 
Recent numerical simulations QMacLow fc Norman 1993| ; |Blondin et al. 1998|) have shown 
that the PDTSO is probably the most dangerous instability, especially for superbubbles 
where the wind shock is adiabatic rather than radiative for a large fraction of its evolution. 
However, the PDTSO works only for the time during which the Mach number of the shell 
is larger than ~ 3. As the PDTSO grows with time as t s with s ~ 1/2, whereas the shock 
velocity decreases as t~ 2 ^ (see equation below) in the radiative phase, amplifications 
30 are typically possible, i.e. the overstability saturates. The above numerical studies 
concluded on this basis that fragmentation of the shell due to these instabilities is unlikely, 
although piondin et al. ( 1998 )| emphasize that in a denser ISM the shock enters the radiative 



phase at a higher Mach number, thus allowing for a longer growth time of the instability. 

Even if the shell can survive dynamical instabilities, it can be fragmented by R-T 
instabilities generated before blowout. As the radiative bubble expands, its pressure drops 
and eventually becomes equal to the that of the ISM. In this case, the wind shock stalls, but 
the shocked wind is still adiabatic and forces the shell to continue its expansion. However, 
the expansion rate is now so slow that the acceleration of the shell is dominated by the 
galactic gravitational field. Thus, the polar region of the shell will be R-T unstable even 
in the absence of a blowout. We calculate the time at which the pressure-confined phase 
occurs (and hence the R-T instability starts to grow) as follows. For constant luminosity, 
the evolution of the shell radius in the radiative phase is given by flCastor, McCray~&| 
I Weaver 19751; IWeaver et al. 19771 Pstriker & McKee 19881) 



-^38 \ ,3/5 



R(t) ~ 66 f 6 /J pc (24) 

where L^s = L/(10 38 erg s" 1 ) is the wind mechanical luminosity, no is the ambient gas 
number density, and te = t/(10 6 yr). From the pressure balance condition R 2 = c 2 ., where 
c s is the ISM effective sound speed, this transition occurs at 



i L \ 1/2 

t p ~ 3 x 10 7 —J- yr (25) 
where c S) io = c s /(10 km s _1 ) is the ISM effective sound speed. The growth time of the R-T 
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instability on spatial scales A s ~ is 

-1.2 

' i I t I 1 / \ i 1 I 'iit! \ 

tRT 



~2ng (A - 1)*" 


-1/2 




[ K (A + i). 




m 



1 Myr (26) 

if the density contrast between the shell and the hot cavity gas A = n s /rih ^> 1 and 
(g) = 10~ 8 cm s -2 . As £rt is much shorter than t p , the shell is fragmented soon after 
it enters the pressure- confined phase. For the case of a Gaussian SFH, the mechanical 
luminosity is initially very low, causing the shock to stagnate much earlier than the 
constant-luminosity case. 

In addition, the shell can be re-accelerated early in its evolutionary stage shortly 
after the first supernovae of the association due to the sudden increase in the mechanical 
luminosity. We have numerically evaluated (g) due to this process assuming the 
Kompaneets approximation is valid, and find that the product of (g) and the duration of 
the re-acceleration (typically ~ 1 — 4 Myr) is ~ 5. Thus, the growth of the instability 
is moderately non-linear. A more detailed study of this process, including the role of 
poisoning, will be given in Dove, Fcrrara, fe Shull (1999) . Therefore, in contrast to the case 



of a constant mechanical luminosity, the shell can be R-T unstable before blowout. In this 
case, a fragmenting shell would allow a higher amount of ionizing radiation to escape the 
H I disk since fragmentation occurs during the time interval in which the photon luminosity 
of the association is near its peak value. 

The shell might also be fragmented by a gravitational instability. Several studies have 
investigated this possibility, both under interstellar (|Elmegreen fc Lada 1977]; fVlcUray ~Sz 



Kafatos 1987| ; [Voit 198*5] ) and intergalactic conditions ( pstriker fc (Jowie 1981] |Couchman 



Rees 1986|) . Here, we simply state the the growth time of the most unstable mode is 



_ 3c lsheii l(] 7 L 1/8 v~ 1/2 vr (27) 

Again, we see that this mechanism is able to fragment the shell on time scales of interest 
and similar to the ones found for the R-T instability. We conclude that, even in the absence 
of blowout, fragmentation of the shell is likely about 3 x 10 7 yr after the beginning of the 
energy injection. 

Even if the shells of superbubbles do not blow out, do they fragment while the 
OB association is still producing ionizing photons? Simple conditions on the occurrence 
of blowout have been derived in different contexts by |Mac Low fc McCray (1988)| and 



Mac Low fc Ferrara (199*9]| . Neglecting the role of magnetic fields, [Koo fc McKee (1992 



studied in detail the case for the Milky Way and concluded that the blowout requires 
L 3S ^, 4. For a marginally confined bubble, we find from equation (|25| ) and equation (]2 
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that t p ~ 6 x 10 7 yr and ta ~ 2.4 x 10 7 yr for the R-T and the gravitational instability, 
respectively. Note that, as L 38 is decreased, fragmentation occurs earlier in time. Thus, 
even bubbles that do not blow out of the disk can fragment on time scales short enough that 
at least some ionizing photons escape. Determining the covering fraction of the fragmented 
clumps is outside the scope of the present estimates and will require numerical simulations 
of the fragmentation process. As a final caveat, we stress that a dynamically important 
magnetic field could qualitatively modify the above conclusions by suppressing some of the 
instabilities (|Tomisaka 1998|) and by introducing new ones. 



Although our models assume that, outside the superbubbles, the neutral gas is 
distributed diffusely, there is observational evidence that a significant fraction of the gas 
is in cold H I clouds, with densities much larger than the space-averaged density. Due 
to the higher recombination rates, these clouds modify our radiation-transfer results and 
therefore merit some consideration. Since most radiation that escapes the H I disk comes 
from young OB associations, with shells having small column densities, we considered the 
ramifications of a clumpy medium by ignoring the effects of evolving superbubbles. A more 
proper treatment would consider the time-dependent evolution of the clouds due to passage 
of the shell of the superbubble as well as the expansion and evaporation of the clouds after 
being photoionized. 

As an approximate model, we let a mass-fraction of the gas, a c \, be contained in cold, 
dense clouds. The remaining fraction of the gas, (1 — a cl ), is distributed diffusely using 
the Dickey-Lockman density distribution (Dickey & Lockman 1990). We consider two 
cloud geometries: (1) spheres of radius r c \] or (2) cylindrical disks with radius r c \ and an 
aspect ratio 5 P = h/r c \, where h < r c \ is the height of the cylinder. Each of the cylindrical 
disks is assumed to be "face-on" with respect to the line of sight from the OB association. 
For either geometry, all clouds are assumed to have identical column densities, sizes, and 
densities. The number of clouds per unit area of the disk is given by 



where V c \ is the volume of each cloud, n c \ is the density of hydrogen within each cloud, and 
Njjj = 6 x 10 20 cm -2 is the observed total column density of hydrogen through the entire 
disk at the solar circle. 



4.2. Clumpiness of the Diffuse H I Gas 



N cl 



a c \N H i 




We assume that the clouds are distributed randomly within the disk of the galaxy 
using a probability distribution, P d (Z), that is constant in the plane and proportional to 
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the Dickey-Lockman density distribution in the vertical direction. The average number of 
clouds intersecting a line of sight having an inclination angle 9 (// = cos 6), emanating from 
an OB association situated on the midplane, is given by 



N _^iV cl _ f ^(S) (ife)" (so^) _1 : for cylindrical disks 

2/1 I °-f(m (^y 1 (so^r ■ ^ spheres 

(29) 

As in §2.3.2, the fraction of LyC photons from a single OB association that escape each 
side of the disk is given by 

tt/2 

^csc(5)= / l^lsm(9)d9, (30) 



where 



Ana 



(2) 



oc 



f(6,S) = l-^- I n 2 H (r)r 2 dr. (31) 



Here, n#(r) is given by the diffuse density distribution for radii outside any cloud and is 
equal to n d for radii inside any cloud. For a given value of 6 and S, we evaluate f(9,S) 
by averaging over 100 Monte-Carlo simulations. For each simulation, the number of clouds 
along the line of sight is drawn from a Poisson distribution having an average number of 
clouds given by equation (|29|) . The distances from the origin of each of these clouds to the 
OB association are picked randomly using the probability distribution P c i(Z). For spherical 
clouds, the average chord length between the line of sight and the cloud is I = h = 4r cl /3. 
For the cylindrical clouds, we assume that each cloud is orientated face-on, such that 
I = 5 p r c \. The observed column density of each cloud in the line of sight is n c \l. Using the 
average values of f(9,S), we determine r] esc (S) via equation (|30|). Finally, as in DS94, we 
determine the total fraction of ionizing radiation that escapes the disk is determined by 
integrating over the luminosity function, 

S JVes C (S)^dS 
(/esc) = 2 . (32) 

J^dS 

Si 

In Table [I], we give (/ csc ) for various values of a c \,r c \, and 8 p . It is apparent that (/ eS c) 
is sensitive to these parameters, which are poorly constrained observationally. For a given 
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Table 1. (/ esc ) for Clumpy Static ISM 



Geometry 


r c \ (pc) 
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<5 P 


S 2 (10 49 s" 1 ) 


(/esc) 


disk 


5 
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value of a c \, (/esc) decreases as the average number of clouds within a line of sight increases. 
We find that the I-front is trapped within most clouds. Only clouds very close to the OB 
association are not able to trap the I-front and become completely ionized, a consequence of 
the high recombination rate within the cloud. Essentially all escaping radiation emanates 
from lines of sight that are free of clouds. For comparison, for the case where a c \ = 0, 
(/esc) = 0.12,0.18,0.22 for S 2 (10 49 s^ 1 ) = 100,250, and 500, respectively (these results 
are identical to those of DS94). We conclude that, if an appreciable amount of the H I 
gas is in cold clouds having a disk geometry and column densities and hydrogen densities 
similar to the H I clouds studied by Fitzpatrick fc Spitzer (1997)| , then (/ esc ) can be a factor 
of ~ 2 — 5 smaller than previously estimated. For a fixed amount of gas locked in cold 
clouds, spherical clouds occupy a smaller solid angle as compared to thin cylindrical clouds, 
allowing more radiation to escape. Note that, since we did not consider superbubbles in 
these calculations, these numbers are only indicative of the relative decrease of (/ eS c) when 
cold clouds are included in the model. 



4.3. Poisoning 

Although the survey used for calibrating the production rate of ionizing photons is 
heavily biased towards associations far away from the parent molecular cloud, most if not 
all associations were born in a molecular cloud. Therefore, early in the lifetime of the 
OB association, it is possible that a considerable amount of matter is photoablated from 
the molecular cloud. This material increases the interior density of the bubble, leading to 
radiative cooling, potentially at a rate equal to the the injected mechanical luminosity, a 
process independently dubbed "poisoning" by [McKee ( 1986 )| and [Shull fc Saken (1995)| . 
This poisoning process can alter the evolution of a young superbubble, as well as the 
propagation of the I-front through the cavity. 

To date, the photoablation rate has not been calculated for a time-varying luminosity 
of an OB cluster with a proper radiation spectrum. In addition, for a given mass-injection 
rate, solving the evolution of the injected gas within the cavity is a complex problem. 
These complexities lead to uncertainties in the radiative cooling rate, both in magnitude 
and in its time evolution. Therefore, in order to estimate the degree to which poisoning 
alters the evolution of the superbubble and (f eS c), we introduced a factor, g, representing 
the fraction of the mechanical luminosity injected into the cavity that is drained due to 
radiative cooling from both photoablated gas and evaporated gas from the dense shell. 
Thus, (1 — g)L mec h is the portion of energy available to do PdV work. In what we consider 



an extreme case, we let g = 0.5 for the duration in which R s h is less than 30 pc [see |Shull 
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& Saken (1995) for a discussion of when cooling due to poisoning becomes negligible]. In 
this case, for our non-coeval model, we find that (/ esc ) increases by roughly 35%. This 
increase occurs because the superbubble shell expands more slowly, and the density of the 
shell is smaller (it is proportional to the Mach number squared of the shell), and because 
the column density of the shell increases more slowly. 

On the other hand, it is possible that enough material is photoablated into the 
hot cavity that, after recombining, it contributes to the absorption of ionizing radiation 
emitted by the OB association, partially offsetting the increase in (/esc) due to the slower 
shell expansion. However, to model this process requires a knowledge of the relative 
configurations of the molecular cloud, the OB association, and the disk of the Galaxy. 
Also, the evolution of the injected gas, especially the expansion of the gas as it expands 
from the molecular cloud or shell boundaries and into the the line-of-sight between the OB 
association and the halo, is a numerical problem beyond the scope of this paper. 



4.4. Triggered Star Formation 

The fraction of ionizing radiation escaping the H I disk would be increased if an 
appreciable number of OB associations were formed within pre-existing dynamic chimneys. 
The two-dimensional filling factor of dynamic chimneys is estimated by 

/mi ~ /asiV" as r ca 7r^ y , (33) 

where / as is the fraction of OB associations that are large enough to produce a dynamic 
chimney, r ca is the lifetime of the cavity, and R cy is the radius of the cavity. As discussed 
in §4.1, associations having an average mechanical luminosity L moc h ^ 4 x 10 38 erg s _1 
are thought to be capable of blowing out of the Galactic disk (|Mac Low fc Ferrara 1999| ; 
|Koo fc McKee 1992| ). This corresponds to an OB association size of iVr ^ 400, or a peak 
luminosity So ^ S CT = 5 x 10 50 s _1 . The fraction of associations that produce a dynamic 
chimney is then 

J dN/dSo (s^ 1 '* 



f ^'t— ^7^7^(1^ r_1 < (34) 



JdN/dS V5l 

Si 

where we assumed S 2 > 5 cr > Si. For T = 2 and Si — 1 x 10 48 s -1 , / as ~ 2 x 10~ 3 . The 
lifetime of the dynamic chimney is estimated to be the sound crossing time, 

-^cy or A 10 km s~ lN \ f (T h 



Tch ~ — ~ 35 ( — —— ) ( p-jjj^; ) M W- (35) 
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Therefore, the probability that an OB association is born within a dynamic chimney is 
given by 

/mi - 0.18 ( e a ( (^-Y, (36) 



6.4 Myr kpc 

and it is unlikely that (/esc) is enhanced significantly by this process if OB associations 
are born randomly in space. However, as suggested by [McCray fc Kafatos 19871 , the birth 
of OB associations could be triggered by the passage of the shock wave associated with 
another superbubble. Indeed, recent far ultraviolet, Ha, and H I observations of several 
dwarf galaxies show evidence of secondary star-formation sites on the dense H I rims of 
superbubbles (|Voit 1988| ; |Stewart 19981) . 



5. Discussion and Conclusions 



We find that (/ esc ) ~ 12 — 20% for the coeval star- formation history and (/ esc ) ~ 4 — 10% 
for the Gaussian star-formation history. For S2 = 2 x 10 51 s _1 , (3 = 1, and T = 2, we 
obtain (/ esc ) = 13% for the coeval SFH and (/ eS c) = 6% for the Gaussian SFH. Adopting 
a production rate [], ^L y c — 4.95 x 10 7 cm -2 s -1 , the corresponding photon fluxes are 
$L y c = 2.1 x 10 6 cm -2 s" 1 for the coeval SFH and $ Ly c = 9.9 x 10 5 cm" 2 s _1 for the 
Gaussian SFH. For comparison, (/ eS c) = 15% using the formalism of DS94, where the 
photon luminosity was assumed constant in time and dynamically evolving superbubbles 
were not considered. In order to determine the degree to which the differences between 
our new results and our old results are due to the treatment of expanding superbubbles 
and dynamic chimneys, we also calculated (/ eS c) for the case where we neglect entirely the 
superbubble dynamics but include the time-varying photon luminosities. For the same 
parameters as given above, we find (/ esc ) = 14% for the Gaussian SFH and (/ eS c) = 16% 
for the coeval SFH, values similar to the old model. These small differences are due solely 
to the new treatment of the luminosity function of associations with time-varying photon 
luminosities, as discussed in §3.3.1. Therefore, the inclusion of the evolution of superbubbles 
reduces the fraction of escaping radiation, as the superbubble shells trap more radiation. 
For coeval SFHs, the inclusion of the dynamics only changes (/ eS c) from 8.2% to 6.5%, since 
most radiation is emitted and escapes before a substantial shell has formed. On the other 
hand, for the Gaussian SFH, (/ eS c) drops from 14% to 6%. This reduction of (/ esc ) occurs 



1 This number is 50% higher than the value used in DS94, due to updated stellar atmosphere models 
whose yield of ionizing photons from OB associations is 50% higher than previously estimated ( Sutherland 
fc Shull 19991 ). 
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because the photon luminosity peaks after a shell has formed, causing fewer photons to 
escape compared to the case when superbubbles are neglected. 

In DS94, we argued that (/esc) must be ~ 14% in order for the Reynolds layer to be 
sustained by radiation from OB associations. However, using COSTAR model atmospheres, 
Sutherland fc Shull ( 199971 found that the yield of ionizing photons from OB associations is 



50% higher than previously estimated, and the production rate of ionizing photons by OB 
associations in the solar vicinity is $L y c — 4.95 x 10 7 cm -2 s^ 1 . Thus, the required value of 
(/esc) to sustain the Reynolds layer is reduced to (/esc) ~ 8%. The flux of escaping radiation 
for the Gaussian SFH is about 25% smaller than this required value. However, as discussed 
in §1, since a portion of the Reynolds layer may reside within the H I disk, it need not be 
sustained solely by the radiation escaping the disk. 

Throughout this paper, we have assumed that the formation rate of OB associations is 
constant in time. For starburst galaxies, the observed star formation rate is much higher 
than the time-averaged value. Therefore, since r] esc (t) is at a maximum shortly after the 
formation of an OB association, starburst galaxies having a morphology similar to our 
Galaxy should have higher values of (/ esc ) than those predicted for our Galaxy. In addition, 
starburst galaxies may have lower H I column densities and lower gravitational fields, also 
contributing to a higher escape fraction of ionizing radiation. These considerations will be 
the focus of our subsequent work. 

If OB associations are indeed the sources of radiation responsible for ionizing the 
Reynolds layer, it is possible that an appreciable amount of radiation can escape the Galaxy 
entirely and contribute to the intergalactic radiation field QBechtold et al. 1987] ; |Giroux 



[Shull 199"7| ). If the average escape fraction from low-redshift, star-forming galaxies exceeds 



~ 5%, then the contribution of ionizing radiation to the intergalactic medium by galaxies 
can exceed the contribution for QSOs ( [Shull 19"9~5| ; |Madau fc IS hull 1996| ; |Giallongo, Fontana 



fc Madau 1997| ). We note, however, that our models predict only the flux of radiation 



penetrating the Galactic halo, before a portion of it is absorbed by the portion of the 
Reynolds layer situated high above the Galactic plane. In order to predict the fraction of 
radiation that escapes the Galaxy entirely, we need a self-consistent model of the density 
distribution of hydrogen gas within the halo, including the contribution due to high velocity 
clouds, whose covering fraction is uncertain. |Murphy, Lockman, fc Savage (1995]| found 



that approximately 37% of 21-cm sightlines near quasars showed HVCs with Nhi > 7x 10 17 
cm -2 . However, observations with the Parkes Multibeam Survey flPutman fc Gibson 1999| ) 
show that the spatial distribution of H I in the HVCs is filamentary, suggesting that Lyman 
continuum could leak through. Furthermore, the actual HVC coverage fraction will vary 
considerably around the Galaxy, depending on proximity to spiral arms and upwelling 
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superbubbles. 

For specificity, if we assume (/esc) — 0.06, then the upward ionizing flux is 
$L y c = 1-5 x 10 6 cm -2 s~ x above the H I disk (Z ^ 1 kpc). If this flux is incident on the 
HVCs, and if we use the relationship between Ha emission and ionizing flux, for optically 
thick clouds and photoionization equilibrium ([Bland-Hawthorn fc Maloney 1999| ), then the 
predicted Ha emission intensity is 

I E = 4.6 niR f in4 $LyC 2 r ) = 0.69 R. (37) 
\10 4 cm~ z s" 1 / 

This intensity is roughly a factor of two larger than the measured intensities from the 
Magellanic Stream ([Wcincr fc Williams 1996" ) and the Smith HVC (piand-Hawthorn et al 



19981) . 



In summary, we have shown that the escape fraction of ionizing radiation from 
star-forming regions of disk galaxies is strongly affected by the presence of superbubbles 
and a cloudy ISM. The shells of the expanding superbubbles quickly trap or attenuate the 
ionizing flux, so that most of the escaping radiation escapes shortly after the formation of 
the superbubble. In our models, the escape fraction above 1 kpc range from (/esc) ~ 3 — 20%. 
For a Gaussian SFH, (/ esc ) is roughly a factor of two lower than the results of Pove fc Shull 



(1994b)| , where superbubbles were not considered. If (/ esc ) ^ 8%, then the upward flux is 
sufficient to sustain the Reynolds layer. We have also examined the radiative transfer in 
a two-phase (cloud/intercloud) medium. We neglect the filling factor of hot, low-density 
gas, consistent with recent O VI studies. If a significant fraction of H I is distributed in 
cold clouds with hh ~ 30 cm -3 and column densities Ah ~ (10 19 — 10 20 ) cm" 2 , (/ esc ) can 
be reduced by a factor of ~ 2 — 5 if the clouds have a disk geometry. Thus, the amount 
of ionizing photons that escape the disk into the IGM is likely to be highly variable, and 
will depend on the luminosity distribution of OB associations, the cloud structure, and the 
massive star formation history. If the escaping LyC radiation is ever observed directly, it 
will confirm these theoretical expectations but also verify the likelihood that stellar ionizing 
radiation plays an important role in reionization of the high-redshift IGM. 
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